load c:\manoj\projects\ace\pred_julia pred_w available_ace_fday missing_data_ace ace_matrix ;
load c:\manoj\projects\ace\Julia_W_new w_climate Julia_W  julia_model_drift fday;    
load c:\manoj\geomag\indices\aplist.mat;
%pred_w = pred_w.*24.366*1e-3; % electric fields in mV/m
w_climate = w_climate.*24.366*1e-3;
%ace_matrix = ace_matrix./15;

% selected_days = [2001,08,13;...
%                 2001,08,21;...
%                 2002,07,29;...
%                 2003,04,05;...
%                 2003,06,10;...
%                 2003,06,21;...
%                 2003,12,11;...
%                 2004,01,31;...
%                 2004,12,16;...
%                 2005,06,04;...
%                 2006,05,12];
            
selected_days =[2001,08,17;...
                2001,10,02;...
                2002,09,13;...
                2003,04,05;...
                2003,06,21;...
                2004,10,30;...
                2005,09,02;...
                2006,01,22;...
                2008,01,12];
            

days_to_plot = datenum(selected_days)-datenum(2000,1,1);

%size of fday is 1002, number of ace overlap is 993 - this is because ace
%does not have data on the last 9 days of julia fdays 
% h=figure(1);
% set(h,'Position',[1 31 1920 1096]);
for kk = 1:length(days_to_plot),
    i = find(floor(fday(:,1)) == days_to_plot(kk));
    L = fday_ap >= fday(i,Julia_W(i).k(1)) & fday_ap <= fday(i,Julia_W(i).k(end));
    mean_ap = mean(ap(L));
    obs_minus_climate = w_climate(i,Julia_W(i).k);
    %obs_minus_climate = obs_minus_climate-nanmean(obs_minus_climate);
    obs_minus_climate_minus_ace = pred_w(i,Julia_W(i).k);
    %obs_minus_climate_minus_ace = obs_minus_climate_minus_ace-nanmean(obs_minus_climate_minus_ace);
    climate_minus_ace = obs_minus_climate-obs_minus_climate_minus_ace;
    ace = ace_matrix(i,:);
    rms_error_1(kk) = sqrt(mean(obs_minus_climate.^2));
    rms_error_2(kk) = sqrt(mean(climate_minus_ace.^2));
    subplot('position',[0.15,0.4,0.8,0.5]);
    time_ax    = fday(i,Julia_W(i).k);     
    plot(time_ax+datenum(2000,1,1),obs_minus_climate,'r','LineWidth',2);
    set(gca,'FontSize',16);
    hold on;
    plot(fday(i,Julia_W(i).k(1):Julia_W(i).k(end))+datenum(2000,1,1),...
    pred_w(i,Julia_W(i).k(1):Julia_W(i).k(end)),'b','linewidth',2);
    date_string = datestr(time_ax(1)+datenum(2000,1,1),29);
    title(sprintf('Date - %s, Ap = %4.1f',date_string,mean_ap));
    datetick;
    axis([datenum(2000,1,1)+Julia_W(i).fday+(8+5)/24 datenum(2000,1,1)+Julia_W(i).fday+(18+5)/24 -inf inf]);
    grid on;
    hold off;
    %xlabel('UT (hours)');
    ylabel('Equatorial EF mV/m');
    set(gca,'XTickLabel',[]);
    subplot('position',[0.15,0.15,0.8,0.2]);        
%   subplot(2,1,2);
    plot(fday(i,:)+datenum(2000,1,1),ace,'k');
   datetick;
    set(gca,'FontSize',16);
    axis([datenum(2000,1,1)+Julia_W(i).fday+(8+5)/24 datenum(2000,1,1)+Julia_W(i).fday+(18+5)/24 -inf inf]);
     xlabel('UT (hours)');
     ylabel('IEF Ey mV/m');
         grid;
         set(gcf,'Position',[ 680   678   824   420])
    pause;
end;
